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1.  INTRODUCTION 


The  objective  of  the  present  investigation  is  to  develop  a  finite  element  model  for  use  in 
contact  problems.  These  problems  occur  when  small  clearance  exists  between  adjacent  parts 
of  the  structure  and  this  clearance  is  closed  during  the  load  cycle.  This  leads  to  contact 
between  materials  which  modified  the  structural  response.  The  structural  contact  problem  is 
essentially  a  nonlinear  phenomenon  as  illustrated  by  the  well-known  Hertz  contact  problem 
(Gaertner  1977).  The  nonlinear  nature  of  the  contact  problem  makes  it  suitable  for  solution  by 
the  finite  element  solution.  The  finite  element  method  has  two  distinct  approaches  for  solving 
contact  problems,  which  include  the  direct  method  and  the  gap  element  method. 

In  the  direct  method,  a  number  of  approaches  can  be  used.  One  such  approach  involves 
formulating  the  problem  with  constraints.  The  solution  is  then  performed  with  a  suitable 
mathematical  programming  technique.  The  direct  method  solves  contact  problems  through 
some  sort  of  special  mathematical  programming  technique.  Hung  and  Saxce  (1980)  describe 
a  method  which  uses  variational  principles  to  solve  contact  problems.  Specifically,  the  contact 
problem  is  equivalent  to  a  minimization  of  potential  energy  with  constraints  on  the 
displacements  so  that  surface  overlap  does  not  take  place.  One  disadvantage  of  this 
approach  is  that  it  is  difficult  to  append  to  existing  finite  element  codes. 

Another  direct  approach  is  to  have  the  finite  element  program  check  for  contact,  and  when 
it  is  established,  the  finite  element  grid  is  then  modified.  The  disadvantage  of  this  approach  is 
that  when  the  contact  occurs,  the  finite  element  grid  changes  and,  consequently,  the  total 
number  of  degrees  of  freedom  is  also  changed.  For  relatively  small  problems,  this  approach 
does  not  lead  to  too  much  difficulty,  but  has  definite  disadvantages  as  the  complexity  of  the 
structure  increases. 

The  second  approach  is  to  use  gap  finite  elements  to  simulate  the  contact  characteristics. 
These  elements  have  to  have  low  stiffness  for  open  gap  and  very  large  stiffness  as  contact 
occurs.  A  number  of  variations  of  this  approach  have  been  tried  (Mahmoud,  Salomon,  and 
Marks  1982;  Stadter  and  Weiss  1979).  One  advantage  of  the  gap  element  approach  is  that 
the  finite  element  grid  does  not  change  as  the  contact  occurs.  The  second  advantage  is  that 
this  approach  can  easily  be  appended  to  existing  finite  element  programs. 
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The  objective  in  the  present  investigation  is  to  incorporate  the  gap  element  model  into  an 
axisymmetric  finite  element  analysis.  The  analysis  into  which  the  contact  capability  is  to  be 
incorporated  is  the  program  SAGA  (Zak  1974,  1975).  SAGA  program  is  an  elastic  code  with 
three  degrees  of  freedom  in  the  displacements.  Consequently,  SAGA  permits  not  only  the 
regular  axisymmetric  solutions,  but  also  torsional  solutions.  The  program  allows  for  quite 
general  loading  conditions  including  concentrated,  distributed,  inertia,  and  thermal  loads. 
Material  properties  are  fully  orthotropic,  allowing  for  axis  of  orthotropy  to  vary  in  the  meridian 
plane  as  well  as  in  helical  direction. 

A  further  objective  is  to  have  the  gap  element  capability  in  an  elastic-plastic  program. 
Consequently,  the  first  step  in  the  investigation  is  to  convert  SAGA  into  elastic-plastic  code. 
The  model  for  an  elastic-plastic  model  has  been  developed  previously  for  finite  elements  in 
connection  with  previous  studies  (Zak,  Craddock,  and  Drysdale  1 979;  Hill  1 948).  The 
conversion  of  SAGA  into  elastic-plastic  versions  involves  two  steps.  The  first  step  is  to 
change  the  program  into  incremental  loading.  The  second  step  is  to  introduce  the  yield 
condition  and  the  calculation  of  incremental  plastic  stress-strain  relations  (Zak,  Craddock,  and 
Drysdale  1979).  In  addition  to  this  capability,  it  is  also  intended  to  have  one  version  of  the 
model  capable  of  elastic-viscoplastic  analysis.  This  modification  is  a  relatively  straightforward 
extension  of  the  elastic-plastic  model. 

The  introduction  of  the  gap  element  capability  into  the  modified  version  of  SAGA  will 
involve  a  number  of  steps.  The  first  of  these  involves  certain  geometric  parameters  which 
describe  the  gap  position  in  the  element  system  as  well  as  the  size  and  orientation  of  the  gap. 
The  position  and  orientation  of  the  gap  will  determine  which  points  can  contact  as  the  gap 
closes  under  load.  The  second  step  is  to  have  the  computer  program  possess  a  systematic 
way  to  check  for  contact.  After  the  contact  has  been  established,  the  program  has  to  modify 
the  stiffness  of  the  gap  elements  to  simulate  closure.  As  two  nodal  points  contact,  they  may 
not  move  relative  to  each  other  perpendicular  to  the  contact  surface.  However,  tangential 
motion  has  to  be  allowed. 

The  present  report  discusses  the  various  steps  necessary  to  achieve  these  objectives. 
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2.  INCREMENTAL  ANALYSIS 


2.1  Load  Increments.  The  development  of  the  incrementally  loaded  finite  element  model 
is  necessary  for  two  reasons.  Even  if  the  material  properties  were  to  be  taken  to  be  linear,  it 
is  necessary  to  increment  the  analysis  for  the  purpose  of  the  contact  solution.  The  contact 
problem,  even  for  elastic  materials,  is  a  nonlinear  phenomenon  and,  therefore,  in  order  to 
perform  this  part  of  the  solution,  the  load  has  to  be  applied  by  increments.  At  each  increment, 
the  contact  conditions  have  to  be  checked,  and  if  new  contact  is  established,  the  cap 
elements  are  modified  for  the  next  increment  of  load. 

Incremental  analysis  is  also  necessary  for  the  pu^ose  of  handling  nonelastic  materials  for 
which  incremental  stress-strain  relations  have  to  he  used.  These  relations  are  load 
dependent,  as  will  be  shown  in  the  next  section. 

2.2  Incremental  Elastic-Plastic  Stress-Strain.  The  present  elastic-plastic  material  will  be 
based  on  Hill’s  criterion  (Hill  1948)  for  orthotropic  materials.  In  addition,  it  is  desired  to  include 
in  the  material  model  the  effects  of  kinematic  hardening.  This  hardening  phenomenon  has 
been  proposed  for  isotropic  materials  (Calladine  1969)  and  extended  to  orthotropic  materials 
(Zak,  Craddock,  and  Drysdale  1979). 

To  describe  the  plastic  work  hardening  behavior  of  a  material,  three  things  are  needed: 

(1)  an  initial  yield  condition 

(2)  a  flow  rule  relating  plastic  strain  increments  to  the  stress  and  the  stress  increment 

(3)  a  hardening  rule. 

For  this  work,  the  Hill’s  yield  criterion  (Hill  1 948)  for  orthotropic  materials  was  chosen. 

This  reduces  to  the  von  Mises-Hencky  yield  condition  if  the  material  is  isotropic.  The  Hill’s 
condition  will  be  discussed  more  fully  later,  but  it  can  be  represented  symbolically  as  follovvs: 

F{a,j)  =  ^  (1) 


at  the  point  where  yielding  occurs. 
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The  Prager  or  the  kinematic  strain  hardening  rule  was  used  (Calladine  1969).  This 
hardening  rule  assumes  that  the  yield  surface  retains  its  initial  size  and  shape  but  undergoes 
a  translation  in  the  direction  of  the  plastic  strain  increment.  After  plastic  flow  begins,  the  yield 
surface  can  be  written  as  the  following: 

F  (atj  -  c t,j)  =  1  ,  (2) 

where  the  a:/  represents  the  translation  of  yield  surface.  The  assumption  that  this  translation 
is  in  the  direction  of  the  plrstic  strain  increment  can  be  written  as  the  following: 


d atj  =  c  dz*  , 


(3) 


where  c  is  a  constant  for  the  material. 

The  flow  rule  chosen  was  also  due  to  von  Mise^;  namely, 

dz*  -  —  d  X  ;  d  X  >  0  . 

3o,y 


(4) 


The  scalar  d  X  i  Equation  4  can  be  determined  by  the  fact  that  the  stresses  remain  on  the 
yield  surface  during  plastic  flow.  This  condition  can  be  expressed  as  the  following: 


5F 


-  ‘S)  -  0 


do(/ 


(5) 


Substituting  from  Equations  3  and  4  into  Equation  5  gives  the  following: 


do  -  c  d  X 

v  .  3o" 


d  F 
do,, 


=  0 


(6) 
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Solving  Equation  6  for  d  X  gives  the  following: 


dX 


1 

-  do, 
c  do,  1 


dF  dF 
dok,  dokl 


(7) 


The  total  strain  increment  is  the  sum  of  the  elastic  and  the  plastic  strain  increments.  Thus, 


dza 

u  '/ 


dt,  -  dza 


(8) 


where  (It,  is  the  total  strain  increment. 


Substituting  for  dz°  in  the  elastic  stress-strain  relationship  gives  the  following: 

do,,  =  EiJU  [dzkl  -  de*p,]  , 


(9) 


where  Eijkl  is  the  elastic  stress-strain  matrix. 


Equation  5  can  be  rewritten  as  the  following: 


« *.')  ¥  - 


do,, 


(10) 


Then,  substituting  Equations  4,  7,  and  9  into  Equation  10  gives  the  following: 


~ifkl 


-  if.  d\ 

OOy, 


-  c  dX 

a°, 


>11-0 

do, 


(11) 


5 


This  can  be  rewritten  as  the  following: 


c  w  3F  c  3F  3F  „  3F  dr  I 

""  "  ^  *  p"  35;  as:  *  35;  5s;  A 


Then,  defining  the  bracketed  term  as  1/D  where  D  is  the  scalar, 


n  P  dF  dF  ^  dF  dF 

dokl  datj  3c,y  d<5,j 


Equation  12  can  be  written  as  the  following: 


dX  =  D  E„w  dek, . 


Substituting  Equation  14  into  Equation  4  gives  a  relation  between  the  total  strain  increment 
and  the  plastic  strain  increment  as  follows: 


dz’  --  D  £„„„  dz. 


<^Ei/  3  ^iikl^zkl  ’ 


where  Cijld  is  defined  from  Equation  1 5. 


Then,  Equation  9  can  be  written  as  the  following: 


3  [  ^ijkl  ~  ^ktmn^mnij  ^Zij  ’ 
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or 


dCij  *  AiJkl  dekl  , 


(18) 


where 


A 


Ijkl 


(19) 


For  the  isotropic  case,  c  can  be  evaluated  from  a  uniaxial  test.  It  can  be  shown  that 


c 


2  E‘Et 

3  E-Et  ' 


(20) 


where  Er  is  the  tangent  modulus,  and  E  is  Young's  Modulus.  In  general,  Erwill  depend  on 
the  stress  state  at  a  given  time.  Thus,  by  knowing  the  elastic  constants  and  the  tangent 
modulus,  the  stresses  can  be  obtained  from  the  total  strains. 


2.3  Application  to  Orthotropic  Materials.  In  the  previous  section,  the  incremental  stress- 
strain  relations  have  been  developed  symbolically  in  terms  of  the  orthotropic  yield  function. 
The  specific  expression  for  the  Hill's  yield  criterion  is  given  by  the  following: 


2  2  2  2  2  2 
r- ,  ,  Oil  °22  °33  °12  °23  °13  (7 

F  (<*//)  =  —  +  —  +  —  +  — T  +  —  +  — T  +  V11  CT22a33 


Y2 

r11 


r  22 


'  33 


Y 2 
*12 


'  23 


Vic 


+  Y22  o,  ,  Ojg  +  Y33  o,  ,  o22 


*  1  , 


(21) 


where  V|i>  v22,  and  V33  are  the  yield  stresses  in  simple  tension  in  the  three  orthotropic 
directions,  and  Y12,  Y23,  and  V13  are  the  corresponding  yield  stresses  in  simple  shear.  The  Y 
terms  are  functions  of  the  yield  stresses.  These  functions  include  the  following: 
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(22) 


Equation  21  can  be  corrected  for  strain  hardening  by  substituting  for  stresses  as  shown  in 
Equations  1  and  2.  Using  Equation  21  after  it  has  been  modified  for  strain  hardening  as 
illustrated  by  Equation  2,  it  is  possible  to  differentiate  the  yield  function  relative  to  stress  (Zak, 
Craddock,  and  Drysdale  1979).  These  relations  then  can  be  used  in  the  previous  section  to 
evaluate  the  incremental  stress-strain  relations  given  by  Equation  18. 


There  remains  the  question  of  evaluating  the  hardening  parameter,  C.  This  was  evaluated 
from  a  uniaxial  test  in  the  isotropic  case,  but  there  are  three  uniaxial  tests  which  can  give 
three  different  results  in  the  orthotropic  case.  Consequently,  in  order  to  evaluate  the 
hardening  parameter  for  orthotropic  materials,  some  limitations  are  placed  on  the  orthotropic 
properties  (Zak,  Craddock,  and  Drysdale  1979).  One  of  these  is  that  the  material  is 
transversely  orthotropic;  that  is,  directions  2  and  3  are  identical.  It  can  be  shown  then  (Zak, 
Craddock,  and  Drysdale  1 979)  that  the  hardening  can  be  related  to  testing  in  the  orthotropic 
direction  1; 


C  = 


2  '  En  T 

3  E„-E11r  • 


(23) 


2.4  Elastic-Plastic  SAGA.  Once  the  incremental  elastic-plastic  relations  have  been 
obtained,  the  next  step  is  to  implement  these  into  the  SAGA  program  to  create  an  elastic- 
plastic  version  of  this  program.  A  flow  chart  for  the  new  version  of  SAGA  is  shown  in 
Figure  1.  This  version  is  labeled  SAGAEP.  The  most  noticeable  difference  between  this 
version  and  the  original  SAGA  is  the  addition  of  two  new  subroutines  ELPLSS  and  YIELD. 
Subroutine  YIELD  is  used  to  determine  orthotropic  yield  conditions,  and  ELPLSS  evaluates 
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incremental  stress-strain  relations  in  the  plastic  region.  Less  obvious  changes  involve 
modification  of  existing  SAGA  subroutines,  and  these  include  MAIN,  QUAD,  STIFF,  and 
STRESS. 

3.  GAP  GEOMETRY 

3.1  Gap  Position.  In  the  formulation  of  gap  elements,  there  are  a  number  of  geometric 
parameters  which  either  have  to  be  set  by  input  or  calculated  intei  _.iy  in  the  program.  The 
first  of  these  parameters  has  to  do  with  defining  the  direction  of  the  gap.  This  direction  is  in 
terms  of  an  l-J  grid.  The  l-J  coordinate  system  will  locate  elements  and  nodal  points.  To 
illustrate  how  nodal  points  are  located,  take  a  typical  element  in  an  l-J  grid  (Figure  2).  This 
element  has  nodal  point  !X(N,1)  located  at  the  comer  that  lies  closest  to  the  origin.  The 
remaining  three  nodal  points  IX(N,2),  IX(N,3),  and  IX(N,4)  are  obtained  by  moving  in  the 
counterclockwise  direction  around  the  element. 

The  direction  of  the  gap  will  depend  on  an  inputted  parameter  called  the  principal  gap 
direction.  As  an  example,  if  the  principal  gap  direction  is  chosen  to  be  1 ,  then  the  gap 
direction  will  be  parallel  to  the  J-axis,  while  the  direction  of  closure  will  be  parallel  to  the  l-axis 
(Figure  3).  If  the  principal  gap  direction  is  set  to  2,  the  direction  of  the  gap  will  be  parallel  to 
the  l-axis,  while  the  gap  closure  will  be  in  the  J-axis  direction  (Figure  4). 

3.2  Determining  Nodal  Pairs.  The  gap  direction  will  determine  which  node  points  will  be 
opposite  one  another  along  the  gap.  This  information  is  important  for  reasons  that  will  be 
shown  in  the  following  sections.  Simply  stated,  these  reasons  include  the  following.  First,  the 
program  must  know  which  node  points  to  pair  and  then  test  for  contact;  and  secondly,  when 
contact  is  established,  the  nodal  point  pairs  will  be  used  to  properly  modify  gap  element 
material  properties.  Referring  to  Figure  3,  if  the  principal  gap  direction  is  set  to  1 ,  then  the 
two  nodal  pairs  for  the  element  are  (IX[N,1],  IX[N,2])  and  (IX[N,3],  IX[N,4]).  Similarly,  for  a 
principal  gap  direction  of  2,  the  two  pairs  are  (IX[N,1],  IX[N,4])  and  (IX[N,2],  IX[N,3]),  as  shown 
in  Figure  4. 

3.3  Gap  Orientation.  Another  geometric  parameter  has  to  do  with  the  orientation  of  the 
gap  elements  relative  to  the  R-Z  axes.  This  orientation  will  be  used  to  transform  nodal 
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positions  and  nodal  displacements  into  a  coordinate  system  where  one  axis  coincides  with  the 
direction  of  the  gap,  while  the  other  axis  coincides  with  the  direction  of  closure.  The  second 
use  of  the  gap  orientation  will  insure  that  gap  material  properties  will  only  stiffen  in  the 
direction  of  closure  once  contact  is  established. 

Once  the  principal  gap  direction  is  determined,  the  angle  of  the  gap  can  be  calculated. 

This  angle  is  defined  in  a  manner  indicated  by  Figure  5.  As  an  example,  take  a  typical 
element  with  a  principal  gap  direction  of  1  and  nodal  points  II,  12,  13,  and  14  [equal  to  IX(N,1), 

IX(N,2),  IX(N,3),  and  IX(N,4),  respectively].  The  angle  p  will  be  calculated  by  first  producing  a 
mean  line  through  the  element  (Figure  6).  The  coordinates  of  points  1  and  2  are, 
respectively,  included  in  the  following: 

R(I1)  +  R(I4)  Z{I1)  *  Z(I4) 

.  ,  _  (24) 

*  (Rmeanl  ,  Zmeanl)  , 


R(I2)  +  R(I3)  Z(I2)  +  Z(I3) 

2  ’  2 

■  ( Rmean2  ,  Zmean2)  . 

Using  Equations  24  and  25,  the  angle  p  can  be  calculated  by  using  the  following: 

P  *  ARCtan  (Rmean2  -  Rmean1 ) 

( Zmean2  -  Zmeanl) 


(25) 


(26) 


This  definition  of  p  will  have  coordinate  axis  N  coincident  with  the  direction  of  closure,  and 
coordinate  axis  T  coincident  with  the  gap  direction. 

If  the  principal  gap  direction  is  equal  to  2,  the  angle  calculation  will  be  the  same  as 
previously  mentioned,  except  that  Equations  24  and  25  will  be  replaced  by  the  following: 

R{I1)  +  R(I2)  Z(I1)  *  Z(I2) 

2  ’  2 

=  ( Rmeanl  ,  Zmeanl)  , 
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(27) 


R{I3)  +  R(I4) 
2 


Z(I3)  *  Z(I4) 
2 


(28) 


*  ( Rmean2  ,  Zmean2)  . 


The  last  geometric  parameter  deals  with  the  distance  between  two  nodal  points  of  a  given 
nodal  pair.  This  distance  in  found  by  using  the  following  relationship: 

gap  size  =  ^/[f?( 1 )  -  R(2)Y  *  [Z(1)  -  Z{2)\ 2  ,  (29) 

where  R(1)  and  Z(1)  are  the  R,  Z  coordinates  of  the  first  node  point,  and  R(2),  Z(2)  are  the 
R,  Z  coordinates  of  the  second  node  point. 

4.  TESTING  FOR  CONTACT  AND  LOAD  FACTORS 

4.1  Testing  for  Contact.  The  gap  finite  elements  are  designed  to  simulate  the  real 
behavior  of  gaps.  Real  gaps  will  have  no  resistance  to  closure  before  opposing  material  faces 
contact.  Consequently,  the  gap  elements  must  possess  very  low  stiffness.  After  surface 
contact  has  been  established,  the  elements  must  now  become  very  stiff  in  order  for  the  two 
surfaces  in  contact  to  have  no  relative  motion  perpendicular  to  the  surface.  The  change  of 
the  gap  element  properties  will  be  initiated  when  the  contact  between  any  nodal  points  is 
achieved.  The  check  for  surface  contact  is  done  by  comparing  the  relative  motion  of  nodal 
pairs  described  previously.  Since  the  analysis  has  already  been  put  in  an  incremental  form 
for  the  elastic-plastic  analysis,  it  becomes  a  simple  matter  to  check  for  contact  at  each  load 
increment. 

Determining  whether  or  not  opposing  node  points  have  contacted  is  done  by  first 
transforming  nodal  displacements  and  nodal  positions  into  a  new  coordinate  system.  This 
coordinate  system  coincides  with  the  gap  direction.  The  coordinate  axis  T  coincides  with  the 
direction  of  the  gap,  and  coordinate  axis  N  coincides  with  the  direction  of  closure.  The 
transformation  of  nodal  displacements  from  the  r,  z,  0  coordinate  system  to  the  N,  T,  S 
coordinate  system  is  of  the  following  form: 
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Uf, 

UT 

•  *  t  T  ]  ' 

uz 

4> 

»  4 

(30) 


where  ur,  ur  and  u9  are  the  nodal  displacements  in  the  r,  z,  and  8  directions,  respectively, 
and  u uv  and  us  are  the  nodal  displacements  in  the  N,  T,  and  S  directions,  respectively. 


The  transformation  matrix  [t]  is  as  follows: 


SIN  p  COS  (5  0 
-COSp  S/A/ [3  0 

0  0  1.0 


(31) 


where  the  angle  (3  is  defined  by  Equation  26.  As  can  be  seen  in  Equation  31 ,  the  coordinate 
axes  0  and  S  will  always  coincide  in  this  transformation. 

The  transformation  of  nodal  position  is  similar  to  the  transformation  of  nodal 
displacements.  This  transformation  is  of  the  following  form: 


(32) 


where  P,  and  Pz  are  the  nodal  positions  in  the  r  and  z  directions  respectively,  and  PN  and  PT 
are  the  nodal  positions  in  the  N  and  T  directions,  respectively. 


The  transformation  matrix  is  as  follows: 


S/A/p  COSp 
-COSP  S/A/p 


(33) 
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The  second  step  in  the  test  for  contact  involves  updating  nodal  positions  in  the  N  direction. 
The  updating  occurs  at  the  end  of  every  load  step;  this  can  be  represented  by  the  following: 

PN  (updated)  =  PN  (original)  +  uN  (total) ,  (34) 

where  PN  (updated)  and  PN  (original)  are  the  updated  and  original  node  positions  in  the  N 
direction,  respectively.  The  term  uN  (total)  is  the  total  displacement  of  the  node  point  in  the  N 
direction  at  the  end  of  the  load  step. 

The  contact  condition  of  an  opposing  node  point  pair  can  be  determined  by  using  the 
updated  nodal  positions.  This  is  done  by  comparing  the  difference  of  the  two  original 
positions  to  the  difference  of  the  two  updated  positions.  This  comparison  can  be  represented 
by  the  following: 

P'n  ( original )  -  P„  ( original ) 

— . - - -  f  (35) 

PN  ( updated )  -  PN  ( updated ) 


where  1  and  2  represent  the  first  and  second  node  points  of  the  node  pair,  respectively.  If 
Equation  35  is  positive,  contact  between  the  opposing  node  points  has  not  taken  place.  If 
Equation  35  is  negative,  however,  the  opposing  node  points  have  overlapped  (i.e.,  crossed 
over  one  another). 

4.2  Load  Factors.  When  real  gaps  close,  the  two  opposing  material  faces  will  never 
overlap  one  another.  For  an  effective  simulation  of  real  gaps,  this  behavior  must  be 
duplicated  in  the  program.  It  is  conceivable  that  a  nodal  pair  that  is  not  contacting  at  the  end 
of  a  given  load  step  can  be  overlapping  at  the  end  of  the  next  load  step.  This  can  occur 
because  the  program  only  checks  for  contact  at  discrete  intervals  and  not  continuously.  The 
program  handles  the  problem  of  overlapping  by  splitting  load  steps.  The  first  step  in  this 
process  takes  all  the  nodal  pairs  that  are  overlapping  at  the  end  of  a  given  load  step  and 
calculates  a  factor  for  each  one  that  depends  on  the  degree  of  overlap. 

The  calculation  of  the  load  factor  is  best  illustrated  by  the  use  of  an  example.  Take  a 
nodal  pair  in  which  the  gap  size  closes  by  three  units  in  load  step  while  the  amount  of 
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overlap  at  the  end  of  load  step  is  one  unit.  The  factor  for  the  previous  nodal  pair  will  be 
calculated  by  the  following: 

1  -  erij£ _  -  .666  -  load  factor  .  (36) 

3  -  gap  closure 


Simply  stated,  if  we  multiply  load  step  load  by  .666,  we  will  have  a  situation  where  the  load  is 
just  large  enough  to  make  the  previous  nodal  pair  contact.  It  is  possible  to  have  a  situation 
where  more  than  one  nodal  pair  is  overlapping;  in  this  case,  the  program  will  determine  the 
smallest  load  factor  and  will  multiply  the  load  by  that  amount. 


When  the  load  step  is  multiplied  by  the  load  factor,  there  remains  a  part  of  the  load  that 
has  yet  to  be  applied.  As  an  example,  if  a  load  factor  of  .6  is  calculated  for  a  load  step  of 
1 ,000  lb,  an  additional  load  of  400  lb  remains.  This  calculation  can  be  represented  by  the 
following  steps: 

1,000  lb  -  original  load  step 
x  .6  -  load  factor 

600  lb  -  reduced  load  step  , 

1,000  lb  -  original  toad  step 

-600  lb  -  reduced  load  step 
400  lb  -  unrun  portion  of  the  load  . 

The  procedure  used  in  the  program  involves  applying  the  residue  force  as  an  additional  load 
step. 

5.  GAP  ELEMENT  STIFFNESS 

5.1  Stiffness  Matrix  for  Gao  Elements.  As  was  stated  previously,  gap  element  material 
properties  will  change  when  contact  occurs.  Before  contact,  the  relative  stiffness  between  the 
opposing  material  faces  is  zero.  After  contact,  the  stiffness  between  the  faces  is  theoretically 
infinite.  For  an  effective  contact  analysis,  this  stiffness  behavior  must  be  simulated  by  the 
program.  With  real  gaps,  the  change  from  a  soft  material  to  a  stiff  material  will  occur  at  all 
areas  of  contact.  Since  the  finite  element  program  only  checks  for  contact  at  discrete 
locations  (namely,  the  location  of  nodal  point  pairs),  this  will  be  simulated  by  a  stiffness 
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change  between  opposing  node  points.  This  simulation  is  good  as  long  as  the  nodal  pairs  are 
not  spaced  too  far  apart. 

In  the  present  model,  gap  material  properties  are  changed  by  modifying  gap  element 
stiffness  matrices  directly.  These  stiffness  matrices  are  calculated  in  a  nodal  system  where 
the  degrees  of  freedom  are  defined  as  in  Figure  7.  The  stiffness  matrices  of  the  three- 
dimensional  quadrilateral  elements  in  the  N,  S,  T  system  have  the  following  form:  S L  where 
i,  j  -  1,12.  As  in  all  stiffness  matrices,  the  rows  of  the  matrix  correspond  to  the  force  in 
question,  and  the  columns  correspond  to  the  degrees  of  freedom.  As  an  example,  force  1  will 
be  related  to  the  12  degrees  of  freedom  by  the  following: 

=  ^1.1  U1  +  ^1,2U2  +  ^1,3^3  +  ^1,4U4 

+  ^1.5U5  +  ^1 ,6a6  +  Si, 7^7 

+  S1i8t/8  +  S,  9  Uq  +  S1 10u10 

+  Sj,t|i/ti  +  S1  it2U12  ,  (37) 

where  the  u  terms  are  the  displacements  corresponding  to  the  indicated  degrees  of  freedom. 

The  stiffness  between  two  degrees  of  freedom  can  be  changed  by  modifying  force 
equations  such  as  Equation  37.  Assume  that  the  relative  stiffness  between  degrees  of 
freedom  1  and  4  needs  to  be  increased.  In  Equation  37,  this  is  done  by  setting  S,  4  equal  to 
the  negative  of  S,  v  and  increasing  the  value  of  S,  r  Equation  37  now  becomes  the 
following: 


^1  =  Si,i(Ui  -  U4)  +  S1  2u2  +  Si  3U3 
+  Si  5U5  +  S1  6u6  +  S,  7u7  +  S1  8u8 

+  S,.9«B  +  S,  ,ou10  +  S, ,nUn  +  S112Ui2  .  (38) 
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In  a  similar  manner,  the  coefficient  S4 ,  is  set  equal  to  the  negative  of  S4  4  in  the  equation  for 
force  4: 

F x  =  S4  4  (u4  - 1/-,)  +  S1  2u2  +  S,  3u3 

+  ^1,5^5  ^1,6U6  +  $1.7^7  +  Sl.BU8 

+  ^1,9^9  +  ^1,10U10  +  ®1.11U11  +  ^1,12^12  •  @9) 

Since  all  stiffness  matrices  are  symmetric,  the  following  is  true: 

S1i4=S4i1.  (40) 


Furthermore,  using  Equation  40, 


(41) 


As  can  be  seen  in  Equations  38  and  39,  the  relative  stiffness  between  the  two  degrees  of 
freedom  is  governed  by  the  magnitude  of  S,  ,  (=S4  4).  The  resuit  of  this  step  is  that  now  the 
relative  motion  of  points  1  and  4  will  be  zero  in  the  normal  direction  A' 


5.2  Transformation  of  Stiffness  Matrices.  One  use  of  gap  element  orientation  involves  the 
transformation  of  stiffness  matrices.  The  program  assembles  the  global  stiffness  matrix  in  the 
r.  z,  and  0  direction  coordinates.  When  gaps  close,  the  change  in  material  properties  occurs 
in  the  direction  of  closure,  which  is  defined  in  the  N,  S,  T  coordinates.  This  problem  is 
handled  by  first  modifying  the  element  matrices  in  the  N,  T,  S  coordinate  system  and  then 
transforming  them  back  into  the  r,  z,  9  coordinate  system. 

The  matrix  transformation  is  given  by  the  following: 


[S'].[H]r[s][H], 


(42) 


where  [s]  and  [S']  are  the  element  stiffness  matrices  in  the  N,  T,  S  and  r,  z,  9  coordinate 
systems,  respectively.  The  transformation  matrix  [w]  is  as  follows: 
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(43) 
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mr 


mr 


mrJ 


where  [r]T  is  the  transpose  of  the  transformation  matrix  defined  in  Equation  31. 


5.3  Gap  Material  Properties  Before  Contact.  An  important  issue  in  this  analysis  involves 
modeling  gap  behavior  before  any  contact  between  opposing  material  faces  takes  place.  In 
the  "real  world,"  gaps  have  no  stiffness  between  opposing  faces  before  contact.  In  the 
program,  this  behavior  is  simulated  by  a  very  soft  gap  material;  the  stiffness  of  the  material  is 
small  but  not  zero.  The  reasons  exist  for  this  nonzero  stiffness.  There  are  a  number  of 
matrix  operations  in  the  program  which  would  fail  if  all  coefficients  of  the  stiffness  matrix  for  an 
element  are  zero. 


5.4  Modeling  Various  Contact  States.  A  gap  element  has  four  different  contact  states  for 
each  value  of  the  principal  gap  direction.  For  an  element  with  a  principal  gap  direction  of  1 
and  opposing  nodal  pairs  (1,2)  and  (3,4)  (shown  in  Figure  8),  the  following  four  contact  states 
exist: 


State  No.  1 
State  No.  2 
State  No.  3 
State  No.  4 


both  node  pairs  not  contacting 

(1.2)  contacting;  (3,4)  not  contacting 

(1.2)  not  contacting;  (3,4)  contacting 
both  nodal  pairs  contacting. 


The  gap  element  stiffness  matrix  will  be  modified  differently  for  each  of  the  four  states. 
Assume  that  the  element  matrix  is  modified  in  the  N,  T,  S  coordinate  system  and  that  the 
three  degrees  of  freedom  at  each  node  are  in  the  N,  T,  and  S  directions,  respectively.  The 
proper  modification  of  the  element  stiffness  matrix  is  as  follows: 


State  No.  1  -  all  matrix  elements  remain  small. 
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State  No.  2  -  set  S, ,  =  S44  =  -S,  4  =  -S4 A 

=  large  number  relative  to  other  matrix  elements. 

State  No.  3  -  set  S77  —  S10 =  -S7  10  =  *Si^7 

=  large  number  relative  to  other  matrix  elements. 

State  No.  4  -  set  ^  —  S44  —  S7  7  =  Siqiq 

=  *^1.4  =  *§4,1  =  '^7,10  =  '^10,7 

=  large  number  relative  to  other  matrix  elements. 

If  the  principal  gap  direction  of  the  element  is  equal  to  2,  the  gap  element  stiffness  matrix 
modification  will  be  similar  to  that  shown  previously.  The  only  difference  between  the  two 
modifications  will  be  that  now  the  two  nodal  pairs  are  (1 ,4)  and  (2,3),  as  shown  in  Figure  9. 
The  four  possible  states  of  contact  are  now  the  following: 

State 
State 
State 
State 

The  proper  modification 
the  following: 

State 

State 

State 

State 

=  '^1,10  =  ‘^10,1  =  *^4,7  =  "S7.4 

=  large  number  relative  to  other  matrix  elements. 


No.  1  -  both  nodal  pairs  not  contacting 
No.  2  -  (1,4)  contacting;  (2,3)  not  contacting 
No.  3  -  (1 ,4)  not  contacting;  (2,3)  contacting 
No.  4  -  both  nodal  pairs  contacting. 

of  the  elements  stiffness  matrix  for  the  four  contact  states  is  now 

No.  1  -  all  elements  in  the  matrix  remain  small. 

No.  2  -  set  S1  ^  =  SW10  =  *S1 10  =  ~S10 1 

=  large  number  relative  to  othe.  matrix  elements. 

No.  3  *  set  S4  4  =  Sj  7  =  -S4  7  —  -S7  4 

=  large  number  relative  to  other  matrix  elements. 

No.  4  -  set  S1  j  =  S44  =  S7  7  =  S10 10 
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5.5  Magnitude  of  Stiffness  Change.  For  the  different  contact  states,  we  increase  the 
magnitude  of  various  matrix  elements.  A  question  remains  of  how  much  to  increase  these 
matrix  elements.  In  theory,  as  long  as  the  on-diagonal  and  off-diagonal  terms  are  equal  and 
opposite,  they  can  approach  infinity  without  affecting  any  other  stiffness  besides  the  relative 
stiffness  between  the  two  degrees  of  freedom  in  question.  In  practice,  however,  as  the 
numbers  approach  infinity,  the  round-off  error  will  cause  a  significant  difference  between  the 
on-diagonal  and  off-diagonal  terms.  The  effect  of  this  error  will  cause  an  artificial 
displacement  constraint  for  the  two  nodal  points  in  the  equation.  By  executing  a  number  of 
test  cases,  it  has  been  found  that  the  change  in  the  magnitude  in  stiffness  should  be  between 
order  of  magnitude  of  2  and  4.  This  order  of  magnitude  applies  relative  to  the  material 
surrounding  the  gap. 

5.6  Transverse  Slip.  The  previous  section  described  the  procedure  developed  for 
modifying  the  stiffness  coefficients  for  the  gap  elements  normal  to  the  gap  direction.  If  only 
these  coefficients  are  modified,  then  the  contact  points  are  relatively  free  to  move  relative  to 
each  other  in  the  direction  parallel  to  the  contact  surface.  However,  in  a  real  situation,  one 
would  expect  that  even  this  motion  would  be  somewhat  restricted  due  to  surface  friction 
effects.  Consequently,  in  the  present  analysis,  a  very  simple  model  is  proposed  to  simulate 
this  effect. 

The  effects  of  friction  will  be  related  to  the  stiffness  in  the  direction  of  closure.  This 
relationship  is  made  through  the  use  of  a  frictional  coefficient.  As  an  example,  if  nodal  pair 
(1,2)  in  an  element  is  contacting,  then  the  following  relationships  will  model  friction: 

^2.2  =  =  '^5,2  =  '^2,5  =  friCt'S,.,  ,  (44) 

^3,3  =  ^6,6  =  '^3,6  =  ’^6,3  =  friCt*Si  ,  ,  (45) 

where  "frict"  is  the  friction  coefficient.  Equation  44  relates  the  stiffness  in  the  direction  of 
closure  to  the  stiffness  in  the  direction  of  the  gap.  Similarly,  Equation  45  relates  the  stiffness 
in  the  direction  of  closure  to  the  stiffness  in  the  tangential  direction. 
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6.  COMPUTER  PROGRAM 


6.1  General  Description.  The  flow  chart  for  the  program  incorporating  the  gap  element 
capability  is  shown  in  Figure  10.  This  program  has  been  designated  as  SAAC  (Structural 
Analysis  Axisymmetric  Contact).  It  is  instructive  to  compare  the  chart  to  that  of  Figure  1  for 
the  elastic-plastic  version  of  SAGA  (SAGAEP).  The  most  obvious  differences  between  SAAC 
and  SAGAEP  are  additions  of  new  subroutines  in  Figure  10.  These  subroutines  include 
GAPANG,  RETRANS,  and  TEST.  The  various  steps  necessary  to  implement  the  gap  element 
capability  will  be  described  in  the  following  sections. 

6.2  Geometrical  Parameters.  As  described  in  the  previous  sections,  the  geometry  of  the 
gap  elements  is  inputted  with  few  parameters,  and  these  are  used  by  the  computer  program 
to  establish  additional  information.  Most  of  the  geometrical  parameters  which  are  established 
in  the  program  are  done  in  the  Subroutine  GAPANG.  This  subroutine  is  called  from  MAIN,  as 
shown  in  Figure  10. 

6.3  Contact  Test.  Testing  for  contact  and  the  calculation  of  load  factors  occurs  in 
Subroutine  TEST.  TEST  is  called  irr  Subroutine  SOLV  after  nodal  displacements  are 
calculated. 

An  assumption  was  made  to  the  load  factor  analysis  to  reduce  execution  time.  The 
assumption  is  that  the  contact  condition  of  nodal  pairs  is  determined  by  the  original  load  step 
and  not  by  the  reduced  load  step.  As  an  example  of  this,  if  a  given  number  of  nodal  pairs  are 
overlapping  at  the  end  of  a  load  step,  the  program  will  then  assume  that  these  nodal  pairs  will 
remain  contacted,  even  though  the  load  step  will  be  reduced.  An  advantage  of  this 
assumption  is  that  original  load  steps  are  only  split  once.  A  disadvantage,  however,  is  that 
the  program  might  predict  gap  closure  at  a  different  load  as  compared  to  a  real  structure  with 
the  same  configuration.  This  problem  can  be  minimized  by  keeping  load  steps  small. 

6.4  Stiffness  Modification.  The  modification  of  stiffness  and  the  friction  analysis  are  both 
done  in  Subroutine  RETRANS.  This  change  of  gap  material  properties  takes  place  in  the 

N,  T,  S  coordinate  system,  and  afterwards,  a  transformation  back  into  the  r,  z,  0  system  is 
done. 
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Subroutine  RETRANS  is  called  in  Subroutine  STIFF  after  the  element  stiffness  matrices  are 
calculated  and  before  the  global  matrix  is  assembled. 

The  Subroutine  RETRANS  is  the  remaining  new  subroutine  to  be  added  to  convert 
SAGAEP  to  SAAC.  In  addition  to  the  new  subroutines,  certain  existing  subroutines  had  to  be 
modified.  The  modified  subroutines  include  MAIN,  QUAD,  SOLV,  and  STIFF. 

7.  ELASTIC-VISCOPLASTIC  OPTION 

7.1  Viscoolastic  Model.  An  additional  objective  of  the  present  investigation  is  to  develop 
an  elastic-viscoplastic  option  for  the  SAAC  program.  Since  this  program  has  already  an 
elastic-plastic  material  model,  it  is  a  relatively  straightforward  step  to  include  viscoplasticity. 
The  model  for  the  elastic-viscoplastic  material  has  been  developed  previously  (Zak,  Craddock, 
and  Drysdale  1979)  in  conjunction  with  one  version  of  the  SANX  program.  A  similar  approach 
will  be  used  here. 

The  elastic-viscoplastic  option  involves  incorporating  the  rate  effect  in  the  stress-strain 
relations.  This  rate  effect  is  incorporated  in  an  incremental  fashion  and  can  be  expressed  by 
the  following  rate-dependent  relation: 


do ij  *  Aijki  dtki  *  tj  de*/*  ,  (46) 

where 

da,"*1  is  an  incremental  change  in  stress  from  time  tn  to  tn+1  ; 
dej/1  is  the  incremental  change  in  total  strain; 

A"jk,  is  the  incremental,  elastic-plastic  material  properties  matrix; 

„Dn.1 

diij  is  the  incremental  change  in  the  viscoplastic  strain  rate;  and 
t|  is  the  material  viscous  coefficient. 

By  using  backward  difference  relations,  it  can  be  shown  (Zak,  Craddock,  and  Drysdale  1979) 
that 
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(47) 


di? 


n*  1 


.  VP 


-  e 


*p 


n- 1 


Then,  Equation  46  becomes  the  following: 


j_n*i  A  n  j  /i»i  _  I  •  vpn 

dO/j  a  Aijk,  dtkl  +  TJ  (B/y  -  tij  ) 


(48) 


By  applying  the  principle  of  virtual  work  to  equation  48  and  making  the  usual  finite  element 
substitutions,  the  following  is  obtained: 


{F}  -  {/"}  -  [*]  {5}  , 


(49) 


where  {fv}  is  the  viscous  force  matrix.  The  {fv}  term  can  be  calculated  for  each  element  as 
the  following: 


{r}an  1[b]t  [{z^n}  -  [ivpn''))  dv  ■  (50) 


and  [/c]  is  the  stiffness  matrix  in  the  following: 


[■ « ]  -  [6]  dV. 


(51) 


Thus,  using  Equations  50  and  51 ,  the  viscous  rate  effects  appear  as  a  modification  to  the 
force  vector  which  can  be  calculated  at  the  start  of  a  load  increment.  This  modification  is 
being  used  in  the  development  of  the  elastic-viscoplastic  option.  In  order  to  convert  elastic- 
plastic  SAAC  program  to  the  elastic-viscoplastic,  changes  to  the  program  were  made  in  MAIN 
and  the  Subroutines  STRESS  and  TRISTF. 

8.  NUMERICAL  EXAMPLES 

8.1  Axisvmmetric  Perfectly  Plastic  Problem.  The  purpose  of  the  first  example  is  to  check 
the  elastic-plastic  capability  of  the  SAAC  program.  This  is  done  by  executing  a  simple 
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example  for  which  an  analytical  elastic-plastic  solution  is  possible  (Shield  and  Ziegler  1958). 
This  example  consists  of  a  uniform,  isotropic  cylinder  subject  to  internal  pressure. 

The  analytical  solution  (Shield  and  Ziegler  1958)  is  valid  for  elastic  perfectly  plastic 
material.  The  finite  element  model,  on  the  other  hand,  is  based  on  incremental  stress-strain 
relation  and  cannot  handle  perfectly  plastic  fesponse.  However,  by  using  very  small  ET,  the 
finite  element  approach  can  simulate  very  closely  perfect  plasticity.  Consequently,  in  this 
example,  the  value  of  ET  was  chosen  to  be  .001 E.  The  analytical  solution  results  predict  the 
plastic  boundary  position  and  the  circumferential  strain  as  a  function  of  internal  pressure. 
Before  discussing  the  comparison  of  these  results  with  the  finite  element  model,  it  may  be 
noted  that  the  analytical  solution  uses  Tresca  yield  criterion. 

Figure  1 1  shows  the  elastic-plastic  boundary  vs.  pressure  curves  for  both  the  numerical 
and  analytical  solutions.  As  can  be  seen,  the  two  solutions  coincide  fairly  well,  except  for  a 
slight  difference  due  to  the  different  yield  criterions.  In  this  problem,  the  Tresca  yield  criterion 
will  predict  yielding  at  a  lower  stress  as  compared  to  the  Von-Mises  yield  criterion  used  in  the 
finite  element  model.  This  fact  will  cause  the  elastic-plastic  boundary  obtained  from  the 
analytical  solution  to  move  outward  at  a  faster  rate. 

The  outer  surface  circumferential  strain  vs.  pressure  curves  for  both  the  numerical  and 
analytical  solutions  are  shown  in  Figure  12.  Again,  the  two  solutions  coincide  closely,  except 
for  a  small  difference  due  to  the  different  yield  criterions.  The  outer  surface  strains  are  higher 
in  the  analytical  solution  because  yielding  occurs  at  lower  loads  as  compared  to  the  numerical 
solution. 

8.2  Flat  Gao — Coarse  Grid  Problem.  The  first  gap  problem  studied  is  an  elastic, 
12-element  configuration  where  the  inner  two  elements  represent  the  gap.  This  configuration 
is  shown  in  Figure  13.  The  configuration  is  loaded  under  external  pressure  such  that 
sometime  during  the  loading  cycle,  node  points  10  and  11  will  come  into  contact. 

Studying  the  displacement  vs.  load  curves  of  nodal  points  10  and  11  reveals  several 
behaviors  (shown  in  Figure  14).  Before  contact,  nodal  point  11  has  a  large  downward 
displacement  per  unit  load.  This  behavior  is  expected  since  the  downward  motion  of  node 
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point  11  is  not  resisted  by  the  gap.  Node  point  10  has  a  relatively  small  upward  displacement 
per  unit  load  before  contact.  This  is  also  expected  since  node  point  10  has  no  direct 
downward  load  to  push  on  it  The  only  loads  affecting  node  point  1 0  are  transmitted  through 
elements  4  and  7  by  shearing  forces. 

After  contact,  the  behavior  of  the  displacement  vs.  load  curves  changes.  Downward 
motion  of  point  1 1  is  now  resisted  by  its  opposing  gap  force.  Node  point  10  is  pushed  down 
by  its  opposing  face  and,  in  turn,  moves  downward  more  rapidly.  As  can  be  seen  in  Figure 
14,  the  relative  motion  between  node  points  10  and  11  is  zero  after  contact.  The  two  node 
points  will  move  together  in  a  manner  that  matches  the  classical  Lame’s  solution  for  thick 
tubes.  This  behavior  is  not  surprising  since  after  contact,  the  configuration  will  act  like  a  solid 
axisymmetric  structure. 

Another  important  thing  to  note  is  that  the  two  node  points  contact  in  the  middle  of  a  load 
increment  This  causes  the  program  to  calculate  a  load  factor  and  then  split  the  load  step. 

The  first  portion  of  the  load  increment  assumes  soft  gap  material  properties,  while  the  second 
portion  assumes  stiff  gap  material  properties. 

Hoop  stress  vs.  load  behavior  was  also  studied  for  this  configuration,  as  shown  in 
Figure  1 5.  Before  contact,  elements  6  and  9  will  carry  a  larger  amount  of  hoop  stress  than 
elements  4  and  7.  This  is  because  the  upper  elements  are  loaded  directly,  while  the  lower 
elements  only  have  loads  transmitted  to  them  by  shearing  faces.  This  behavior  is  in  direct 
opposition  to  Lame’s  solution,  wherein  the  lower  elements  have  the  highest  hoop  stresses. 

After  contact,  the  hoop  stress  vs.  load  curves  will  be  similar  to  the  curves  for  a  solid 
structure.  As  can  be  seen  in  Figure  15,  the  hoop  stress  vs.  load  curves  after  contact  have  the 
same  slope  as  their  corresponding  Lame’s  solution  curves. 

8.3  Flat  Gap— Fine  Grid  Problem.  A  configuration  with  the  same  dimensions  and  loading 
condition  as  the  previous  example  was  analyzed.  The  only  difference  between  the  two 
examples  is  the  refinement  of  the  element  grid.  This  new  configuration  is  shown  in  Figure  16. 
Using  this  example,  the  effect  of  the  grid  on  the  contact  analysis  can  be  studied. 
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The  load  vs.  displacement  curves  for  nodal  points  58  and  59  are  shown  in  Figure  1 7. 
These  two  nodal  points  will  contact  at  higher  load  than  their  corresponding  nodal  points  in  the 
coarse  grid  configuration.  This  is  probably  due  to  the  additional  elements  which  cause  the 
finite  element  analysis  to  model  pressure  distribution  over  the  structure  more  realistically. 

The  hoop  stress  vs.  load  curves  is  shown  in  Figure  18.  The  curves  exhibit  the  same 
behavior  as  in  the  coarse  grid  structure,  except  for  a  small  difference  due  to  the  refinement  of 
the  mesh. 

8.4  Angled  Gap  Problem.  This  example  is  included  to  illustrate  the  capability  of  handling 
gaps  which  are  not  necessarily  oriented  along  any  particular  axis.  This  example  involves  an 
elastic,  12-element  configuration  where  the  gap  is  angled  with  respect  to  the  Z-axis.  The 
configuration  and  type  of  loading  are  shown  in  Figure  19.  The  results  are  shown  in  Figure  20, 
which  shows  the  gap  size  vs.  load  curve  at  the  location  of  the  nodal  pair  (10,1 1).  As  can  be 
seen,  the  closure  rate  of  the  gap  is  linear  with  pressure  until  contact  is  established.  After 
contact,  the  gap  size  will  remain  zero  for  the  remaining  portion  of  the  loading  cycle.  This 
behavior  is  what  is  expected  from  this  sort  of  configuration  under  the  indicated  loading. 
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Figure  1.  Flow  Chart  of  SAGAEP. 
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Figure  4.  Gao  Element  Configuration  for  a  Gao  in  I  Director 


Figure  5.  Definition  of  the  Gap  Element  Orientation  Angle. 
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Figure  12.  Outer  Circumferential  Strain  as  a  Function  of  Pressure  for  the  Perfectly  Plastic 
Disc  (E  =  Elastic  Modulus^  “ 
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Denotes  Nodal  Number 


Figure  13.  Numerical  Example  With  Flat  Gap  and  Coarse  Grid. 


PRESSURE,  PS I  t  103 


15.  Hood  Stress  as  a  Function  of  Pressure  for  the  Flat  Gap  Example  (Coarse  Grid) 
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Figure  19.  Numerical  Example  With  Angled  Gap. 
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\PPENDIX: 

INPUT  CARDS  FOR  GAP  C  :PUTER  PROGRAMS 
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1.  TITLE  CARD 


Format  (20A4) 

Columns  1-80 

2.  CONTROL  CARD 

Format  (615,  F5.0,  615) 
Columns  1-5 

6-10 

11-15 

16-20 

21-25 

26-30 

31-35 

36-40 


41-45 

46-50 

51-55 

56-60 

61-65 


TITLE  (Title  for  particular  case) 


NNLA  (Number  of  nonlinear  approximations; 

NNLA  =  1  for  this  version  of  the  program.) 

NUMTC  (Number  of  temperature  cards;  if  -2, 
a  constant  temperature  is  specified.) 

NUMMAT  (Number  of  different  materials;  6  maximum) 

NUMPC  (Number  of  boundary  pressure  cards; 

200  maximum) 

NUMSC  (Number  of  boundary  shear  cards; 

200  maximum) 

NUMST  (Number  of  boundary  shear  cards  in 
tangential  direction;  200  maximum) 

TREF  (Reference  temperature) 

INERT  (This  parameter  decides  if  inertia  loads  will  be 
present;  INERT  +  0  means  zero  values  of  axial 
acceleration  and  angular  acceleration  and  velocity  for 
each  load  increment.) 

NLINC  (Number  of  load  increments  with  time; 

NUNC  ^  1) 

INC  I  (If  INCI  =  0,  then  inertia  loads  for  each  time 
increment  will  be  equal  to  the  inertia  load  for  the  first 
time  increment) 

INCF  (If  INCF  =  0,  then  surface  loads  for  each  time 
increment  will  be  the  same  as  for  first  increment.) 

IPLOT  (Plot  parameter;  1  if  plot  required) 

LGAP  (If  LGAP  =  1 ,  then  the  gap  direction  will  be 
parallel  to  the  J-axis;  if  LGAP  =  2,  then  the  gap 
direction  will  be  parallel  to  the  l-axis.) 
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2a.  TIME  INCREMENT  CARD 


Note:  This  card  is  used  only  for  elastic-viscoplastic  option — omit  for  elastic-plastic 
version. 

Format  (El  0.0) 

Columns  1-10  DELTIM  (Time  increment) 

3.  MESH  GENERATION  CONTROL  CARD 

Format  (515) 

Columns  1-5  MAXI  (Maximum  value  of  I  in  mesh;  25  maximum) 

6- 1 0  MAXJ  (Maximum  value  of  J  in  mesh;  1 00  maximum) 

11-15  NS  EG  (Number  of  line  segment  cards) 

1 6-20  NBC  (Number  of  boundary  condition  cards) 

21-25  NMTL  (Number  of  material  block  cards) 

4.  LINE  SEGMENT  CARDS 

The  order  of  line  segment  cards  is  immaterial,  except  when  plots  are  requested;  in  this 
case,  the  line  segment  cards  must  define  the  perimeter  of  the  solid  continuously.  The 
order  of  line  segment  cards  defining  internal  straight  lines  is  always  irrelevant. 

Format  [3(213,  2F8.3),  15] 

Columns  1-3  I  coordinate  of  1st  point 

4-6  J  coordinate  of  1st  point 

7- 14  R  coordinate  of  1st  point 

t 

15-22  Z  coordinate  of  1st  point 
23-25  I  coordinate  of  2nd  point 
26-28  J  coordinate  of  2nd  point 
29-36  R  coordinate  of  2nd  point 
37-44  Z  coordinate  of  2nd  point 
45-47  I  coordinate  of  3rd  point 
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■  48-50  J  coordinate  of  3rd  point 

51-58  R  coordinate  of  3rd  point 

59-66  Z  coordinate  of  3rd  point 

67-71  Line  segment  type  parameter 

If  the  number  in  column  71  is 

0  Point  (input  only  1st  point) 

1  Straight  line  (input  only  1st  and  2nd  points) 

2  Straight  line  as  an  internal  diagonal  (input  only  1st 
and  2nd  points) 

3  Circular  arc  specified  by  1st  and  3rd  points  at  the 
ends  of  the  arc  and  2nd  point  at  the  midpoint  of  the 
arc 

4  Circular  arc  specified  by  1st  and  2nd  points  at  the 
ends  of  the  arc  with  the  coordinates  of  the  center  of 
the  arc  given  as  the  3rd  point  (delete  I  and  J  for  3rd 
point) 

5  Straight  line  as  boundary  diagonal  for  which  I  of  1st 
point  is  minimum  for  its  row  and/or  I  of  2nd  point  is 
minimum  for  its  row  (input  only  1st  and  2nd  points) 

6  Straight  line  as  boundary  diagonal  for  which  I  of  1st 
point  and/or  2nd  point  is  maximum  for  its  row  (input 
only  1st  and  2nd  points) 

NOTE:  In  specifying  a  circular  arc,  the  points  are  ordered  such  that  a  counterclockwise 
direction  about  the  center  is  obtained  upon  moving  along  the  boundary. 

5.  BOUNDARY  CONDITION  CARDS 

Each  card  assigns  a  particular  boundary  condition  to  a  block  of  elements  bounded  by  II, 
12,  J1,  and  J2.  For  a  line,  II  »  12,  or  J1  =  J2.  For  a  point,  II  =  12,  and  J1  =  J2. 

Format  (415, 110,  3F10.0) 


Columns 

1-5 

Minimum  1 

6-10 

Maximum  1 

11-15 

Minimum  J 

49 


16-20  Maximum  J 


21-30  Boundary  condition  code 

31-40  Radial  boundary  condition  code,  XR 

41-50  Axial  boundary  condition,  XZ 

51-60  Tangential  boundary  condition,  XT 

If  the  number  in  Columns  21-30  is 

0  XR  is  the  specified  R-load, 

XZ  is  the  specified  Z-load,  and 
XT  is  the  specified  d-load. 

1  XR  is  the  specified  R-displacement, 

XZ  is  the  specified  Z-load,  and 

XT  is  the  specified  6-load. 

2  XR  is  the  specified  R-load, 

XZ  is  the  specified  Z-displacement,  and 
XT  is  the  specified  6-load. 

3  XR  is  the  specified  R-displacement, 

XZ  is  the  specified  Z-displacement,  and 
XT  is  the  specified  6-load. 

4  XR  is  the  specified  R-load, 

XZ  is  the  specified  Z-load,  and 
XT  is  the  specified  6-displacement. 

5  XR  is  the  specified  R-displacement, 

XZ  is  the  specified  Z-load,  and 

XT  is  the  specified  6-displacement. 

6  XR  is  the  specified  R-load, 

XZ  is  the  specified  Z-displacement,  and 
XT  is  the  specified  6-displacement. 

7  XR  is  the  specified  R-displacement, 

XZ  is  the  specified  Z-displacement,  and 
XT  is  the  specified  6-displacement. 

6.  MATERIAL  BLOCK  ASSIGNMENT  CARD 

Each  card  assigns  a  material  definition  number  to  a  block  of  elements  defined  by  the  I,  J 
coordinates. 

Format  (515,  2F  10.0,’  215) 
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Columns  1-5  Material  definition  number  (1  through  6) 

6-10  Minimum  I 

11-15  Maximum! 

16-20  Minimum  J 

21-25  Maximum  J 

26-35  Material  principal  property  inclination  angle  BETA 

which  defines  N-S  plane  orientation  relative  to  z 
direction  (see  Figure  4) 

36-45  Material  principal  property  inclination  angle  ALPHA 

which  defines  the  orientation  of  the  N-T  plane  relative 
to  the  r-z  plane  (see  Figure  4) 

46-50  IANG  (If  IANG  »  0,  then  ALPHA  is  same  for  total 

material  block,  if  IANG  -  1,  the  ALPHA  varies  in  sign 
in  the  I  direction  from  element  to  element  every  NANG 
.  elements.  This  will  allow  for  equal  but  opposite  helical 
angles.) 

51-55  NANG  (Number  of  elements  in  the  I  direction  with  the 
same  ALPHA) 

NOTE:  Gap  elements  will  be  identified  by  a  material  definition  number  of  2. 

7.  PLOT  TITLE  CARD 

Format  (20A4) 

Columns  1-60  Title  (Title  printed  under  each  plot) 

8.  PLOT  GENERATION  INFORMATION  CARD 

Format  (2F10.0) 

Columns  1-10  RMAX  (Maximum  r  coordinate  of  mesh) 

1 1-20  ZMAX  (Maximum  z  coordinate  of  mesh) 

NOTE:  Use  only  if  IPLOT  »  1  (plot  required). 

9.  TEMPERATURE  FIELD  INFORMATION  CARDS 

If  NUMTC  In  columns  6-10  of  the  CONTROL  CARD  is  greater  than  1,  the  temperature 
field  is  given  on  cards.  One  card  must  be  supplied  for  each  point  for  which  a 
temperature  is  specified. 
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Format  (3F10.0) 

Columns  1-10  R  coordinate 


1 1-20  Z  coordinate 
21-30  Temperature 

If  NUMTC  in  columns  6-10  of  the  CONTROL  CARD  is  -2,  a  constant  temperature  field 
is  specified;  the  value  is  given  on  a  single  card. 

Format  (FI  0.0) 

Columns  1-10  Temperature 

10.  MATERIAL  PROPERTY  INFORMATION  CARDS 

Cards  10  to  1 1  (1  la  for  elastic-viscoplastic  option)  must  be  specified  for  each  material 
(maximum  of  6). 

a.  MATERIAL  IDENTIFICATION  CARD 
Format  (215,  2F10.0) 

Columns  1-5  Material  identification  number 

6-1 0  Number  of  temperatures  for  which  properties  are 

given  (12  maximum) 

1 1-20  Mass  density  of  material  (if  required) 

21-30  Thermal  expansion  parameter  (if  1,  free  thermal 

expansion  on  the  material  property  cards;  otherwise, 
coefficients  of  thermal  expansion  are  on  the  material 
property  cards.) 

b.  MATERIAL  PROPERTY  CARDS 

Two  cards  are  required  for  each  temperature. 

First  Card 
Format  (7F10.0) 

Columns  1-10  Temperature 

1 1-20  Modulus  of  elasticity,  EN 
21-30  Modulus  of  elasticity.  Eg 
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31-40 

Modulus  of  elasticity,  ET 

41-50 

Poisson’s  ratio, 

51-60 

Poisson’s  ratio,  v^ 

61-70 

Poisson’s  ratio,  vST 

Second  Card 

Format  (6F10.0) 

Columns  1-10 

Shear  Modulus, 

11-20 

Shear  Modulus,  GST 

21-30 

Shear  Modulus,  G^ 

31-40 

«nT  ora* 

41-50 

OgT  ora*. 

51-60 

Oj-T  or  aj- 

YIELD  STRESS  CARDS 

(not  needed  in  elastic  version) 

Format  (7F10.0) 

Columns  1-10 

Yield  stress  in  tension  in  N  direction 

11-20 

Yield  stress  in  tension  in  S  direction 

21-30 

Yield  stress  in  tension  in  T  direction 

31-40 

Yield  stress  in  shear  in  NS  direction 

41-50 

Yield  stress  in  shear  in  NT  direction 

51-60 

Yield  stress  in  shear  in  TS  direction 

61-70 

Hardening  parameter;  C 

11a.  VISCOSITY  CARD 

NOTE:  This  card  used  only  for  elastic-viscoplastic  option;  omit  for  elastic-plastic 
version. 
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Format  (El  0.0) 

Columns  1-10  ETA  (viscosity  coefficient) 

INERTIA  LOAD  CARD 
Format  (3F10.0) 

Starting  with  this  input  card  and  including  the  boundary  force  cards,  this  data  is  to  be 
input  as  a  block  for  each  load  step,  that  is,  NUNC  times.  There  are  the  following 
exceptions  to  this: 


Columns 


(a)  If  INERT  =  0,  then  this  card  is  to  be  omitted 
completely  (no  inertia  load). 

(b)  if  INCI  =  0,  then  this  card  is  not  repeated,  but  appears 
in  first  block  only  (the  inertia  loads  are  constant  for 
each  load  step). 

(c)  If  INCF  =  0,  then  the  following  boundary  pressure  and 
shear  cards  are  to  be  given  only  for  the  first  block  and 
not  repeated  again  (the  pressure  and  shear  loads  are 
constant  for  each  load  increment). 

1  -1 0  ACELZ  (axial  acceleration) 

1 1-20  ANGVEL  (angular  velocity) 

21-30  ANGACC  (angular  acceleration) 


13.  BOUNDARY  PRESSURE  CARDS 

One  card  is  required  for  each  boundary  element  which  is  subjected  to  a  normal 
pressure;  that  is,  the  number  of  these  cards  is  NUMPC  for  each  load  increment. 

Format  (215,  FI 0.0) 

Columns  1-5  Nodal  point  M 

6-1 0  Nodal  point  N 

11-20  Normal  pressure 

As  shown  in  Figure  A-1,  the  boundary  element  must  be  on  the  left  when  progressing 
from  M  to  N.  Surface  normal  tension  is  input  as  a  negative  pressure. 


14 *  SOUNDARV  gHEAR  naano 

•ha,  is,  ■»*» — n 

Format  (215,  FI  0.0) 

Columns  f_=  ...  .  . 

'-o  Nodal  point  M 

6-10  Nodal  point  N 

11-20  Surface  shear 

ss:  t£sr.tar  »  -  ^ 


Fi9Ure  A'2-  Example  of  Boundary  Shear 


Loadinr 
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15.  BOUNDARY  TRANSVERSE  SHEAR  CARDS 

One  card  is  required  for  each  boundary  element  which  is  subject  to  transverse  shear; 
that  is,  the  number  of  these  cards  is  NUMSC  for  each  load  increment. 


Format  (215,  FI  0.0) 

Columns  1-5  Nodal  point  M 

6-1 0  Nodal  point  N 

1 1-20  Surface  transverse  shear 


Figure  A-3.  Example  of  Boundary  Transverse  Shear  Loading. 
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interest  for  which  the  report  will  be  used.)  _ _ — - 
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